CO2 fertilization of terrestrial photosynthesis inferred from site to global scales

Significance The magnitude of the CO2 fertilization effect on terrestrial photosynthesis is uncertain because it is not directly observed and is subject to confounding effects of climatic variability. We apply three well-established eco-evolutionary optimality theories of gas exchange and photosynthesis, constraining the main processes of CO2 fertilization using measurable variables. Using this framework, we provide robust observationally inferred evidence that a strong CO2 fertilization effect is detectable in globally distributed eddy covariance networks. Applying our method to upscale photosynthesis globally, we find that the magnitude of the CO2 fertilization effect is comparable to its in situ counterpart but highlight the potential for substantial underestimation of this effect in tropical forests for many reflectance-based satellite photosynthesis products.

V egetation photosynthesis is responsible for the largest flux of carbon from the atmosphere into the biosphere (1). Both theory and experimental observations show that the resulting gross primary production (GPP) increases with rising atmospheric CO 2 , a process known as the CO 2 fertilization effect (CFE) (2)(3)(4)(5)(6)(7)(8). The CFE plays a strong role in offsetting anthropogenic emissions by directly increasing the terrestrial carbon assimilation rate. It affects both the global carbon (6) and water budgets (9), which leads to changes in temperature and precipitation, and indirectly affects climate change through vegetation-climate feedbacks (10) by increasing plant growth (11). Understanding the CFE is thus critical to understanding the evolution of Earth's climate.
Despite the importance of global photosynthesis, however, there are few long-term records of observationally inferred GPP in natural environments. This leads to large uncertainty in the magnitude of photosynthetic change over time and the CFE (4,(12)(13)(14). Eddy covariance (EC), a measurement of gas exchange deployed at hundreds of ecosystems worldwide, provides estimates of photosynthesis that could contain information on the CFE (14). Estimating CFE from such observations is challenging, however, because of the large confounding effects from climatic variability of natural environments and the short duration of measurements at many sites. Studies that have attempted to isolate the CFE from carbon fluxes have been limited either to simple statistical approaches that preclude causative attribution (15,16) or to process-based approaches that, despite capturing the sign of responses, report widely varying sensitivities (17,18). This uncertainty has led to debate regarding the magnitude of the CFE, as evidenced by the large spread between satellite-and Earth system model (ESM)-based CFE estimates (19)(20)(21)(22)(23). Indeed, because of the lack of observational evidence at the ecosystem scale, many GPP products based on satellitederived metrics do not explicitly account for CFE (20,24,25), while current representations of CFE also vary markedly among ESMs (26)(27)(28).
Here, we leverage globally distributed EC observations and reanalysis data to detect and attribute CFE using an eco-evolutionary optimality (EEO) framework. This framework reconciles three wellestablished optimization theories in combination with Fick's law and the Farquhar-von Caemmerer-Berry (FvCB) photosynthesis model (29), constraining photosynthesis and key intermediate variables, namely stomatal conductance (g) (30,31), intercellular leaf CO 2 concentration (c i ) (32,33), and photosynthetic capacity (34)(35)(36)(37)(38) (see Materials and Methods for general descriptions and SI Appendix for derivations). The framework can analytically diagnose GPP and its sensitivity to seven measurable variables, namely atmospheric CO 2 concentration (c a ), leaf area index (LAI), air temperature (T a ), volumetric soil water content (SWC), specific humidity (q a ), incident shortwave radiation (SW in ), and surface pressure (P), without the need to prescribe biome-specific photosynthetic properties (35). Our results show an overall increase in photosynthesis at measurement sites worldwide and attribute a large proportion of the increase to the effect of rising CO 2 .

Significance
The magnitude of the CO 2 fertilization effect on terrestrial photosynthesis is uncertain because it is not directly observed and is subject to confounding effects of climatic variability. We apply three well-established eco-evolutionary optimality theories of gas exchange and photosynthesis, constraining the main processes of CO 2 fertilization using measurable variables. Using this framework, we provide robust observationally inferred evidence that a strong CO 2 fertilization effect is detectable in globally distributed eddy covariance networks. Applying our method to upscale photosynthesis globally, we find that the magnitude of the CO 2 fertilization effect is comparable to its in situ counterpart but highlight the potential for substantial underestimation of this effect in tropical forests for many reflectance-based satellite photosynthesis products.
Strong CO 2 -Induced Increase in Site-Scale Photosynthesis. Three metrics are used to quantify and assess the response of GPP to CO 2 .
• β CO2 (gC m À2 year À1 ppmv À1 ) is the partial differential sensitivity of GPP to CO 2 for given environmental conditions (i.e., ∂GPP ∂ca in Eq. 1). It is directly diagnosed by the EEO framework, which decouples the confounding effects on GPP from other environmental conditions. • β ln app is a dimensionless, apparent logarithmic GPP response ratio to CO 2 following Ref. 8 for comparison with other studies (Eq. 3). Changes in GPP include both CO 2 and non-CO 2 effects.
• β ln dir is a dimensionless, direct logarithmic GPP response ratio to CO 2 , which means changes in GPP include CO 2 effects only. Confounding effects are neglected or analytically decoupled through ideal experimental manipulations such as free-air CO 2 enrichment (FACE) experiments or partial derivative approaches. It shares the same equation as β ln app , and is only used for comparison with other studies.
We find a strong increasing trend of GPP measured across 632 site-years of observations in the EC network of 9.  (41). Since β ln app is derived from natural environments and influenced by CO 2 and climate variability, it differs from the direct β ln dir discussed later. The detected overall trends are robust to the uncertainty caused by filtering of low-quality data, and the uneven distribution of sites and site-years (SI Appendix, Fig. S2). Variations in the EC-inferred GPP due to partitioning methods and friction velocity filtering methods have a minimal effect on the EEOinferred GPP (shaded area in Fig. 1A), where the former is used to calibrate the canopy upscaling of the latter. In general, the EEO-inferred GPP reproduces the interannual variability (IAV) of the EC-inferred GPP (Pearson's r = 0.91) ( Fig. 1A; SI Appendix, Fig. S1). The high level of agreement from annual ( Fig. 1B) to monthly (SI Appendix, Fig. S3) scales supports the use of the EEO framework for attribution analysis.
To quantify the GPP trends contributed by different factors, we use the EEO framework to perform a univariate sensitivity analysis (SI Appendix, Fig. S4). The aggregated trend from individual factors (10.3 gC m À2 year À2 , IQR ∈ [7.7, 11.0]) is similar to the EEO-inferred and EC-inferred trends ( Fig. 2A). We find a strong direct contribution from c a at the site level, which accounts for 44% (4.5 gC m À2 year À2 , P < 0.001) of the aggregated GPP trend ( Fig. 2A). We also diagnosed another overall CO 2 -induced GPP trend using the partial differential approach (4.9 gC m À2 year À2 , the first term on the right-hand side of Eq. 1), resulting in a small difference (0.4 gC m À2 year À2 ) compared with the univariate analysis. We then convert the estimate from the univariate analysis, which corresponds to a direct response ratio (β ln dir ) of 0.61 for GPP to CO 2 . This β ln dir is lower than those light-saturated leaf-level analyses in FACE experiments (β ln dir = 0.79) (4,8) and those using deuterium isotopomers (β ln dir = 1.0) (8,42) but is comparable to a theoretical β ln dir of 0.6 in Ref. 8 that considers canopy radiative transfer, suggesting a coexisting of light-saturated and light-limited photosynthesis over our study sites. CFE is stronger under light-saturated conditions but is present regardless of light conditions (27,28). An increase in c a eventually elevates leaf c i through optimization (SI Appendix, Eq. S16), which in the FvCB scheme leads to higher light use efficiency and carbon assimilation rates (7,32).
In addition to c a , T a plays the second most important role, accounting for 28% (2.9 gC m À2 year À2 , P < 0.05) of the aggregated trend ( Fig. 2A), mainly because of the effect of warming on photosynthetic capacity, i.e., the maximum rate of Rubisco activity, V cmax , and the maximum rate of electron transport, J max (Fig. 2C) (38). Although increasing temperature negatively affects GPP through, for example, vapor pressure deficit (VPD) and some biochemical reaction pathways (Fig. 2C), these negative effects are weaker than the positive effects of temperature ( Fig. 2 A and B). SWC and specific humidity, as proxies for plant water supply and atmospheric water demand, together contribute 14% of the aggregated trend (0.91 and 0.50 gC m À2 year À2 , P > 0.05) but are statistically insignificant. These two factors (i.e., SWC and q a ), however, are the largest contributors to GPP IAV, and their contributions (median values are 48 and 33 gC m À2 year À1 , respectively; Fig. 2B) are an order of magnitude larger than those of the GPP trends ( Fig. 2A). These results are consistent with previous data and ESM-driven studies on plant responses to water (43)(44)(45)(46)(47) and are further supported by observations that GPP anomalies are regulated by stomatal conductance and leaf c i , as plants attempt to maximize their efficiency of carbon gain per unit water loss (30,32,48,49). Other factors, such as long-term changes in LAI and radiation, could also lead to changes in GPP across the EC network. However, we find a small and statistically insignificant trend in GPP due to increased LAI (0.88 gC m À2 year À2 ; Fig. 2A), which is consistent with a previous study that solely considered the response of GPP to changes in LAI (50). This small LAIinduced trend in GPP is comparable to the LAI trend, both of which are on the order of & year À1 (51), although our framework diagnoses a large underlying sensitivity of GPP to LAI (mean ∂GPP ∂LAI = 119, intersite SD = 178, unit: gC m À2 year À1 ). In addition, light saturation in canopy radiative transfer (52) may lead to a low LAI-induced GPP trend as leaf area increases do not effectively translate to increases in the fraction of absorbed photosynthetically active radiation at high LAI (3,20). For the changes in incoming shortwave radiation, their contribution to the aggregated GPP is also small (0.43 gC m À2 year À2 ; Fig. 2A) since there is no large trend in radiation over the study sites at the annual scale. Therefore, we conclude that neither changes in LAI nor radiation significantly contribute to the observed upward trend in GPP. The effect of changes in surface pressure is also negligible ( Fig. 2 A and B).

A B
Effects of CO 2 Fertilization at the Global Scale. To examine global CFE patterns, we apply our framework to Moderate Resolution Imaging Spectroradiometer (MODIS) LAI and European Centre for Medium-Range Weather Forecasts Reanalysis Version 5 -Land data (hereafter, ERA5-land). Here, our EEO framework is calibrated by the ensemble mean of all eight satellite-derived GPP products. We then assess both the GPP trend directly estimated by the EEO framework and the underlying differential sensitivities, which themselves can be used to reconstruct a composite trend estimate (Eq. 1). The two approaches show almost identical apparent GPP trends of about 4.7% decade À1 (P < 0.001) relative to their mean GPP throughout 2001 and 2016 (A0 and A1 in Fig. 3A). The resulting CO 2 -induced GPP trend is about 4.1% decade À1 (4.4 gC m À2 year À2 , P < 0.001) (A2 in Fig. 3A; SI Appendix, Fig. S5). We note that the fractional contribution of c a to the apparent GPP trend at the global-scale analysis may be uncertain, and we do not report it because of the uncertainties in apparent GPP trend caused by the uncertain trends in climate forcings. However, estimates of the CO 2 -induced GPP trends (A2 in Fig.  3, calculated as the product of β CO2 and Δc a ) are independent of the uncertainties in the climate forcing trends because of the nature of our partial differential approach: 1) The spatial variation of β CO2 is rational because the climate forcings well capture their spatial variations at the global scale (53-55), and 2) Δc a is derived from the National Oceanic and Atmospheric Administration (NOAA)'s observations with high confidence. As a result, after excluding the spatial heterogeneity introduced by LAI and canopy structure, we show a good agreement when comparing the leaf-level CO 2 -induced GPP trend at the globalscale analysis with their EC site counterparts (r = 0.88; SI Appendix, Fig. S6A).
We compare our CO 2 -induced GPP trends (A2 in Fig. 3) with simulations from 13 dynamic global vegetation models (DGVMs) and eight satellite-derived GPP products during their overlap period (2001 to 2016). Our estimates are at the upper end of their trend distribution, at least one-third stronger than their median (Fig. 3A). The lower GPP trends in the median of DGVMs (B1) and satellite-derived products (C1) is likely due to their smaller CFE in evergreen broadleaf forests (EBFs) (Fig. 3B). Our EEO framework suggests similar relative GPP trends caused by CO 2 among all C3 species (including EBFs), ranging from 4.3 to 5.1% decade À1 (A2 in Fig. 3B). We further replace the ensemble mean of all satellite-derived GPP with each individual product to calibrate the EEO framework, finding a small spread in the diagnosed CO 2 -induced relative  Fig. S3D). Attribution of GPP trends is computed by univariate sensitivity analysis. The other variables are kept as the climatological mean of their study period. Blue * and × indicate statistically significant (P < 0.05) and insignificant trends in the estimated GPP, respectively. The hollow circles, thick black vertical lines, and gray horizontal lines are the median, IQR, and mean of trends. "EC-inferred": GPP trends estimated from the EC network. "EEO-inferred": GPP trends directly estimated by the EEO framework. "Aggregated": sum of GPP trends for each factor from the sensitivity analysis. (B) GPP IAV is proxied by the product of partial derivatives of GPP to control factors and IAVs of control factors. Partial derivatives are diagnosed using climatological means of the control factors and IAVs are defined as 1 SD of the data with long-term trends removed. (C) GPP IAV induced by temperature IAV through different pathways. K is the Michaelis-Menten coefficient, Γ* is the CO 2 compensation point, and η* is the ratio of water viscosity under ambient air temperature to a reference temperature at 298.15 K. In B and C, each violin consists of values from 68 sites.
GPP trends (SI Appendix, Table S2), despite the fact that these individual satellite GPP products used for calibration vary markedly (global GPP in 2001 with a minimum of MODIS Collection 5.5 [MOD-C55] at 104 Pg C year À1 and a maximum of Pmodel-s0 at 131 Pg C year À1 ). In EBFs, it is noteworthy that the satellite GPP products that only consider the indirect CFE on satellite reflectance, but not the direct CFE on gas exchange [i.e., FluxCom, MOD-C55, MOD-C6, and vegetation photosynthesis model (VPM)], do not report an increasing GPP trend (Fig. 3B). Their GPP trends are not a proxy of CFE because they model GPP solely based on climate forcings and leaf greenness or its surrogates. The climate forcings in reanalysis data may show negative impacts on GPP because of excessive temperature and drought stress in the tropics (see the EBF column in SI Appendix, Table S3), and the weak trends in leaf greenness cannot fully extrapolate the direct CFE effect because of the saturation of reflectance in EBFs (56)(57)(58)(59)(60)(61)(62). While the choice of climate forcing could introduce some uncertainties to our CFE estimates, we find similar diagnostic results with another forcing dataset [i.e., Climatic Research Unit and Japanese 55-year Reanalysis (CRU-JRA55) in SI Appendix, Fig. S7]. Furthermore, we calculate β ln app and β ln dir for the globalscale analysis (SI Appendix, Table S3), which results in the same conclusions as those using the relative GPP trend metric.
The Partial Differential Sensitivity of GPP to CO 2 . To further characterize the underlying CFE, we show the spatial distribution of β CO2 , which disentangles the confounding effects from other factors. The magnitude of β CO2 is largest in hot and humid regions at the annual scale (Fig. 4A). Since the CO 2 trend is assumed to be constant worldwide, the spatial pattern of β CO2 represents that of CO 2 -induced GPP trends (A2 in SI Appendix, Fig. S5). At the ecosystem level, β CO2 is highest in EBFs at about 5.8 gC m À2 year À1 per ppmv CO 2 , followed by croplands and other forests from the temperate to boreal regions at about 2.2 gC m À2 year À1 per ppmv CO 2 (Fig. 4B). EBFs have the longest growing season, and their high temperature and radiation determine the high photosynthetic capacity (V cmax and J max ) and thus β CO2 (4,28,35). On the other hand, croplands cultivated in mild regions also have high photosynthetic capacity (63) and additional resource inputs from human management, including multiple cropping and irrigation (51), leading to a relatively high annual β CO2 . Compared with C3 plants, our results suggest a positive but much weaker β CO2 in C4 plants (3,4,13), on the order of 0.6 gC m À2 year À1 per ppmv CO 2 , because the C4 CO 2 pump inhibits the benefit of CFE (49).
According to plant photosynthesis optimization theories and gauged by the GPP magnitude estimated from an EC network or satellites, our EEO framework analytically constrains the sensitivity of GPP to CO 2 (i.e., β CO2 ) from the gas exchange perspective. Our results suggest that the analytical solution of β CO2 is insensitive to the GPP uncertainty used for calibration of the framework (Fig. 1A; SI Appendix, Table S2) and largely unaffected by interannual fluctuations of climatic conditions in recent years (Fig. 4B). This highlights the advantage of the EEO framework in that the partial derivatives decouple the confounding effects of climate on CFE. A recent study using regression methods showed a 40% decline in CFE from 1982 to 2015 (16). However, their key findings are subject to the large temporal uncertainty in the satellite data and limited by the causal ambiguity of their regression (64)(65)(66). In addition, our EEO framework assumes that plants coordinate nutrient investment to optimize photosynthetic capacity (V cmax and J max ) with environmental constraints (34-37), rather than prescribing . "A0": estimated directly from our EEO framework. "A1": diagnostic trends composited from partial derivatives using the EEO framework. "A2": diagnostic CO 2 -induced trends through the partial derivative approach using the EEO framework. "B1": 13 DGVM models, all forcing time-varying. "B2": 13 DGVM models, CO 2 only. "C1": eight satellite-derived GPP products. For B1, B2, and C1, the bars represent their median. The gray bars represent GPP trends considering all effects. The red bars represent GPP trends caused by CO 2 only. BEPS, Boreal Ecosystem Productivity Simulator; BESS, Breathing Earth System Simulator. specific values for plant functional types. Photosynthetic capacity not only controls the slope of the GPP response to CO 2 but also strongly influences it by determining whether photosynthesis is light saturated or light limited (27,28). We explore the influence of coordination at different timescales, comparing the results with the EC-inferred GPP ( Fig. 1A; SI Appendix, Fig.  S8), and adopt a decadal coordination to the environments in the peak LAI month for nutrient investment (67, 68) while adjusting the enzyme kinetics monthly with temperature (38). These assumptions result in a fixed reference J max /V cmax ratio for each site (or grid) over the study period while still allowing both to vary spatially to account for different growing climates. The decadal coordination allows for the presence of lightsaturated and light-limited photosynthesis because the coordination is set at a much longer timescale than the stomatal and leaf c i optimization (35). Under these conditions, the EEO framework diagnoses a large proportion of marginally lightsaturated photosynthesis, which responds more strongly to changes in c a than does light-limited photosynthesis (27,28). Specifically, 63 to 69% of C3 canopy photosynthesis originates from light-saturated conditions for all seasons using the ERA5land forcing. Together, these features explain the high β CO2 predicted by our framework and also lead to a good agreement with the overall trend in EC-inferred GPP and supported by multiple independent studies (39)(40)(41)(42).

Conclusion
We detect a large increase in GPP from a globally distributed EC network. Our EEO framework successfully captures the trends and IAVs in the EC-inferred GPP. The site-scale analysis suggests that CO 2 is the dominant contributor to the GPP trends, while SWC and specific humidity control the IAVs. Using this framework to scale the GPP globally, the estimated absolute CO 2 -caused GPP trend is comparable to its EC-inferred counterpart and translates this CO 2 fertilization effect to a global increase in photosynthesis of 4.1% decade À1 since the 2000s (or 5.1 PgC decade À2 for global sum, 4.4 gC m À2 year À2 for global average). However, our EEO framework diagnoses a high CFE compared with the median values of the DGVMs and the satellite-derived GPP products, with the largest disagreement in tropical forests. These tropical forests are responsible for about one-third of global GPP, and thus, an accurate estimation of CFE in tropical forests is critical to global carbon cycle modeling. Finally, we urge expansion of the currently sparse EC observing network in tropical ecosystems to improve understanding of their physiological responses to CO 2 and climate change.

Materials and Methods
EEO Framework for GPP and CFE. We derive a parsimonious EEO framework that constrains plant photosynthesis and water loss at a monthly scale. In addition to attributing changes in GPP through factorial simulations, the framework allows for partial differential sensitivity of GPP to atmospheric CO 2 and various measurable environmental conditions. This is particularly important as neither GPP nor CFE can be measured directly in the natural environment, and the apparent changes in GPP are due to the confounding effect of changes in CO 2 and other environmental conditions. We use Fick's law and FvCB photosynthesis model as the governing equations for leaf-level photosynthesis. Photosynthesis is then constrained by three well-established optimization theories, which are widely supported by theoretical and empirical evidence (3, 4, 30, 32, 34-37, 44, 48, 69-73). Reconciling these optimization theories aims to constrain key variables in photosynthesis (33), such as the marginal water-use efficiency (mWUE) (74), stomatal conductance (g), leaf intercellular CO 2 concentrations (c i ), and photosynthetic capacity. Details of the EEO framework are derived in the SI Appendix, Text S1 to S7. We briefly summarize the concept of this framework below.
According to Fick's law of mass transfer, neglecting the leaf boundary layer and mesophyll resistances, leaf-level photosynthesis is the product of stomatal conductance (g) and the gradient between atmospheric (c a ) and leaf intercellular CO 2 concentrations (c i ). For a first intuitive guess, an increase in c a can lead to increased photosynthesis. However, plants actively optimize the tradeoff between carbon gain and water loss. Optimization partially offsets the direct photosynthetic benefits of elevated c a by indirectly reducing g and increasing c i to increase the carbon gain per unit water loss eventually.
• Theory 1 states that vegetation maximizes the sum of carbon gain and water loss through optimizing g in response to changing environment (30). • Theory 2 states that vegetation minimizes the summed cost of transpiration and carboxylation per unit carbon assimilation through optimizing c i (32).
Specifically, theory 1 constrains g, and theory 2 constrains c i . Combining theories 1 and 2 provides a direct constraint on the mWUE. The optimization of g and c i should be sufficient and necessary because not only can the variation in g affect c i by controlling the inflow rate of ambient CO 2 (30), but also the variation in c i can regulate stomatal aperture by adjusting guard cell membrane potential (4,75,76). To further consider the water stress on photosynthesis, we use a logistic function to approximate the change of mWUE with SWC (SI Appendix, Fig. S9) (33,49,72,77). This directly allows the framework to incorporate the impact of SWC into the photosynthetic optimization processes (SI Appendix, Fig. S10) by calibrating a site-specific constant (i.e., ζ o ; SI Appendix, Text S4).
In terms of the biogeochemical response, we assume that plants optimize nutrient allocation in response to the elevated CO 2 and their growing environment. We use the photosynthetic coordination theory (theory 3) to constrain how plants optimize their nutrient investment, which allows for the estimation of photosynthetic capacity (i.e., V cmax and J max ) without prescribing biome types (34)(35)(36)(37). bisphosphate carboxylase (Rubisco) activity and ribulose-1,5-bisphosphate (RuBP) regeneration (interchangeable with light saturated and light limited, respectively) (34)(35)(36)(37).
In other words, the coordination theory states that leaf nutrient content is regulated to reflect photosynthetic capacity, which is opposite to the hypothesis that photosynthetic capacity is limited by leaf nutrient content (78). Although the theory is well supported by multiple studies (34)(35)(36)(37), determining the causation between the photosynthetic capacity and leaf nutrient content, or indeed a third alternative of finding common factors that limit the two, warrants further investigation (78). Nevertheless, the coordination theory captures the correlation between photosynthetic capacity and leaf nutrient content and implies that nutrient limitation should be implicit in LAI (37). We tested different timescales of the photosynthetic capacity coordination with respect to the monthly optimization of g and c i and compared our results with the EC-derived GPP (SI Appendix, Fig. S8). We use the timescale that best captures the trend and IAV of the EC-derived GPP (Fig. 1A) as the basis of our analysis. That is, the reference rates of V cmax and J max acclimate to the average environment of the peak LAI month during the study period (67,68). For the other months, V cmax and J max are a function of temperature and their reference values (38). Therefore, our framework allows the photosynthetic capacity to be optimized to the changing atmospheric CO 2 and meteorological conditions on a decadal timescale, and the photosynthesis is roughly colimited by Rubisco and RuBP. Other timescales explored include reference V cmax and J max values optimized according to 1) year-toyear variation of peak month CO 2 + average of peak month meteorological conditions during the study period, 2) year-to-year variation of peak month CO 2 and meteorological conditions, 3) month-to-month variation of CO 2 + average of peak month meteorological conditions during the study period, and 4) month-to-month variation of CO 2 and meteorological conditions. These represent a range of acclimation timescales, from shortterm monthly and annual acclimation to longer-term decadal acclimation.
In order to preserve the parsimonious and analytical nature of the EEO framework, leaf-level photosynthesis is upscaled by LAI using a big-leaf approach. Changes in LAI implicitly account for canopy-scale effects of growing season extensions, recovery from disturbance, and succession. The upscaling factor for light extinction, denoted as j G μ ð Þ μ j, represents the monthly average canopy shape (G) and solar zenith angle (μ) of a particular location (SI Appendix, Text S5). This factor is obtained by calibration with the EC-inferred GPP for site-scale analysis and satellite-derived GPP for global-scale analysis. j G μ ð Þ μ j for each month is then kept constant without year-to-year variation. We stress that the statistical significance of annual trends and interannual variations in EEO-inferred GPP are not caused by calibration but by the seven forcing variables: ambient CO 2 concentration (c a ), LAI, air temperature (T a ), volumetric SWC, specific humidity (q a ), incident shortwave radiation (SW in ), and surface pressure (P). We finally evaluate the EEO framework qualitatively, which meets six of the seven criteria proposed by Ref. 72 (SI Appendix, Text S6).
A key benefit of the analytical framework, in addition to allowing a direct estimate of GPP, is that it can be used to assess the partial derivatives of GPP to different variables, which can further be used to reconstruct a composite estimate of GPP changes for attribution analysis. The following equation can express the change in GPP and sensitivities of GPP to multiple factors: where ΔGPP is the GPP change; ∂GPP ∂ca Δc a is the direct CFE; ∂GPP ∂ca is the partial derivative of GPP to CO 2 , which reflects the absolute response of GPP to unit change in CO 2 ; and so on. We treat these variables as independent forcings to different physiological processes within the EEO framework, but these physiological processes eventually interact to influence the gas exchange. For those larger-scale land-atmosphere interactions, the covariances between the forcing variables are implicit in their measurements (or reanalysis data). We denote ∂GPP ∂ca as β CO2 for simplicity. The symbol Δ can either be the interannual variation or long-term trend. To decompose the sensitivity of GPP to T a into different pathways (i.e., Fig. 2C), we additionally use the chain rule of calculus, as follows.
where D is VPD, V cmax should be replaced with J max if the photosynthesis is light limited, K is the Michaelis-Menten coefficient, Γ* is the CO 2 compensation point, and η* is the ratio of water viscosity under ambient air temperature to a reference temperature at 298.15 K. The full analytical form of the partial derivatives can be retrieved from the MATLAB script in the SI Appendix.
To compare with other studies, we compute the changes in GPP relativized by changes in CO 2 with a metric used in Ref. 8: where GPP s and c a,s are of values in the starting year and GPP e and c a,e are of values in the ending year. A β ln value of a unity indicates direct proportionality GPP response to CO 2 . We remind that this β ln factor is only for intercomparison purposes; it is not a partial derivative and does not, by default, isolate the response of GPP to factors other than CO 2 as β CO2 does. Since our framework can analytically attribute the trends to individual factors, we distinguish direct CO 2 responses, β ln dir , from apparent "CO 2 responses," β ln app .
FLUXNET Dataset. In this study, site-scale EC-inferred GPP and meteorological forcings are provided by the FLUXNET2015 Tier 1 dataset, and the uncertainties associated with this dataset have been previously reviewed (14). To calculate the EEO-inferred GPP (and its sensitivities to various factors), we use air temperature "TA_F," incoming shortwave radiation "SW_IN_F," surface pressure "PA_F," VPD "VPD_F" (convertible to specific humidity), and volumetric SWC "SWC_F_MDS_1" in the dataset as inputs. The EEO framework performs calculations on a monthly basis for daytime averages.
1. Filtering of meteorological data. We only use data in which air temperature, incoming shortwave radiation, and VPD are measured or gap filled with good quality [i.e., quality flag (QF) is 0 or 1]. The same filtering is also applied to the EC-inferred GPP (i.e., FLUXNET2015 GPP). 2. Calculation of monthly daytime averages. We first compute the monthly average of the variables at a given half-hour/hour and aggregate them to monthly daytime averages. Daytime is determined by the nighttime flag (QF = 0) provided by the dataset and further checked with the half-hourly/ hourly shortwave radiation (>0 Wm À2 ). 3. Calculation of reference EC-inferred GPP. In order to compare our results with EC-inferred GPP, we use the mean of "GPP_NT_VUT_MEAN," "GPP_DT_VUT_MEAN," "GPP_NT_CUT_MEAN," and "GPP_DT_CUT_MEAN" as the main reference of EC-inferred GPP. In addition, we also use the individual EC-inferred GPP based on different partitioning and friction velocity filtering methods to analyze the corresponding GPP uncertainty propagation (shaded area in Fig. 1A). 4. Exclusion of low-quality monthly data. We exclude all the monthly averages with less than 50% good-quality timestamps or negative GPP. 5. Aggregation of annual GPP. Since the quality filtering from steps 1 to 4 may lead to data gaps, in annual aggregations, we ensure that no sites have more than one missing value of EEO-inferred GPP and EC-inferred GPP during the growing season (defined as monthly climatological mean EC-inferred GPP > 30 gC m À2 mo À1 ). If there is one and only one such missing value during a year, EEO-inferred and EC-inferred GPP will be filled by their corresponding low-quality data without the filtering of step 4. Otherwise, that particular year will be excluded. 6. Calculation of GPP annual anomalies. We compute the annual anomalies of the EEO-inferred GPP and EC-inferred GPP for each site independently across the observational network (SI Appendix, Fig. S1) and aggregate the anomalies across sites (Fig. 1A). 7. Last, our evaluations are performed for those sites with more than 5 y of good-quality data at the annual scale.
After the preprocessing, there are 68 qualified flux sites (632 site-years) from 2001 to 2014, covering 10 different biome types, including croplands (10), closed shrublands (1), deciduous broadleaf forests (11), EBFs (4), evergreen needleleaf forests (18), grasslands (12), mixed forest (6), open shrublands (2), savannas (1), and woody savannas (3) (SI Appendix, Table S1). In the site-scale analysis, all biome types are treated as C3 species because of lack of information. The GPP trends of site-scale analysis are examined by the Mann-Kendall test (R package: https://cran.r-project.org/web/packages/zyp/index. html). An overall uncertainty of the site-scale GPP trends is characterized by the IQR of trends through resampling of site-years ( Fig. 2A; SI Appendix, Fig. S3D). In addition, we calculated the IQRs due to site heterogeneity by resampling of sites (SI Appendix, Fig. S3 B and E) and due to the uneven distribution of sites and site-years over the study period by randomly shuffling the first-year data (SI Appendix, Fig. S3 C and F). The shaded area in Fig. 1A shows the variations in the EC-inferred GPP due to the partitioning and frictional velocity filtering methods and their propagation to the EEO-inferred GPP. On the other hand, we analyze the dataset without the quality filtering in steps 1, 2, and 4. This allows a comparison of our EEO-inferred GPP with the fully gap-filled FLUXNET data (SI Appendix, Fig. S3A), ensuring that the resulting findings are not due to data filtering.
We estimate another EEO-inferred GPP using CRU-JRA55 forcing (2001 to 2018, 0.5°× 0.5°, https://catalogue.ceda.ac.uk/uuid/13f3635174794bb98cf8ac 4b0ee8f4ed) for comparison. We use the following variables: "temperature at 2m," "maximum temperature at 2m," "downward solar radiation flux," "pressure," and "specific humidity." Since there is no volumetric SWC in CRU-JRA55, we adopt it from ERA5-land. Daytime temperature is averaged from "temperature at 2m" and "maximum temperature at 2m." MODIS Leaf Area Index. We use Collection 6 MODIS LAI as our LAI inputs (52). We refine the MODIS Terra and Aqua LAI (8-d frequency, 500 m) by checking their quality flags following a previous study (51). The quality of this dataset has been extensively validated (55) and reported elsewhere (51,79). We convert the original 8-d LAI into monthly frequencies. For the flux-site analysis, we match the 500-m LAI to each site's longitude and latitude and average over in a 3 × 3 window (1.5 km × 1.5 km). For the global-scale analysis, we convert the 500-m LAI into 0.5°× 0.5°.
Atmospheric CO 2. We use the monthly average atmospheric CO 2 concentration from the Mauna Loa Observatory and the South Pole Observatory provided by NOAA's Earth System Research Laboratory (https://www.esrl.noaa. gov/gmd/ccgg/trends/). DGVM Gross Primary Production. DGVM outputs from the Trends in Net Land-Atmosphere Exchange (TRENDY-v6) project are used in this study (http://dgvm. ceh.ac.uk/node/9) (1). Here, we use two model experiments from TRENDY-v6, either varying CO 2 only (time-invariant "preindustrial" climate and land use mask, S1) or varying CO 2 , climate, and land use (all forcing time-varying, S3). We exclude Sheffield DGVM (SDGVM) as it has an unrealistic sudden drop in GPP after 2007. More details are documented in SI Appendix, Table S4.  Table S4.
Land Cover Maps. In order to analyze the spatial pattern of our diagnosed results, we use Collection 6 MODIS yearly product as the reference land cover map (2001 to 2018, 0.05°× 0.05°) (85). We refine the yearly maps for the study period into a single map and resample to 0.5°× 0.5°by taking the mode class of each grid cell. We aggregate the International Geosphere-Biosphere Program classification types into five C3 biomes (EBFs, other forests, short woody vegetation, grasslands, croplands) and one C4 biome (SI Appendix, Fig. S11). There is no biome aggregation for EBFs. Other forests are aggregated from deciduous broadleaf forests, mixed forests, evergreen needleleaf forests, deciduous needleleaf forests, and woody savannas. Short woody vegetation includes closed shrublands, open shrublands, savannas, and permanent wetlands. Croplands include croplands and cropland/natural vegetation mosaics. Finally, grids with greater than 50% C4 vegetation are labeled as C4 according to the ISLSCP II C4 vegetation percentage (https://daac.ornl.gov/cgi-bin/ dsviewer.pl?ds_id=932).